import numpy as np
import pandas as pd
import geopandas as gpd
import rasterio
from rasterio import features
import matplotlib.pyplot as plt
import sklearn
from sklearn.metrics import confusion_matrix, f1_score, accuracy_score, classification_report
from pathlib import Path
from IPython.display import display
import plotly.express as px
import plotly.figure_factory as ff
import plotly.graph_objects as go
import plotly.offline
plotly.offline.init_notebook_mode()
print('All libraries successfully imported!')
print(f'Scikit-learn: {sklearn.__version__}')
All libraries successfully imported! Scikit-learn: 0.24.2
computer_path = '/export/miro/ndeffense/LBRAT2104/'
grp_letter = 'X'
lut_path = f'{computer_path}data/LUT/'
# Directory for all work files
work_path = f'{computer_path}GROUP_{grp_letter}/WORK/'
in_situ_path = f'{work_path}IN_SITU/'
classif_path = f'{work_path}CLASSIF/'
am_path = f'{work_path}ACCURACY_METRICS/'
Path(am_path).mkdir(parents=True, exist_ok=True)
print(f'Accuracy Metrics path is set to : {am_path}')
Accuracy Metrics path is set to : /export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/ACCURACY_METRICS/
site = 'NAMUR'
year = '2020'
feat_nb = 2
no_data = 0
reclassif_flag = True
if reclassif_flag:
field_classif_code = 'grp_A_nb'
field_classif_name = 'grp_A'
else:
field_classif_code = 'grp_1_nb'
field_classif_name = 'grp_1'
s4s_lut_xlsx = f'{lut_path}crop_dictionary_new.xlsx'
in_situ_val_shp = f'{in_situ_path}{site}_{year}_IN_SITU_ROI_VAL.shp'
in_situ_val_tif = f'{in_situ_path}{site}_{year}_IN_SITU_ROI_VAL.tif'
if reclassif_flag:
classif_tif = f'{classif_path}{site}_{year}_classif_RF_feat_{feat_nb}_{field_classif_name}.tif'
else:
classif_tif = f'{classif_path}{site}_{year}_classif_RF_feat_{feat_nb}.tif'
cm_csv = f'{am_path}{site}_{year}_CM.csv'
cm_html = f'/export/miro/ndeffense/LBRAT2104/GIT/eo-toolbox/notebooks/7_Classification/figures/{site}_{year}_CM.html'
am_html = f'/export/miro/ndeffense/LBRAT2104/GIT/eo-toolbox/notebooks/7_Classification/figures/{site}_{year}_AM.html'
print(f'Classification file used : {classif_tif}')
print(f'Validation polygons used : {in_situ_val_shp}')
Classification file used : /export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/CLASSIF/NAMUR_2020_classif_RF_feat_2_grp_A.tif Validation polygons used : /export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/IN_SITU/NAMUR_2020_IN_SITU_ROI_VAL.shp
print(f'Raster template file : {classif_tif}')
# Open the shapefile with GeoPandas
in_situ_gdf = gpd.read_file(in_situ_val_shp)
# Open the raster file you want to use as a template for rasterize
src = rasterio.open(classif_tif, "r")
# Update metadata
out_meta = src.meta
out_meta.update(nodata=no_data)
crs_shp = str(in_situ_gdf.crs).split(":",1)[1]
crs_tif = str(src.crs).split(":",1)[1]
print(f'The CRS of in situ data is : {crs_shp}')
print(f'The CRS of raster template is : {crs_tif}')
if crs_shp == crs_tif:
print("CRS are the same")
print(f'Rasterize starts : {in_situ_val_shp}')
# Burn the features into the raster and write it out
dst = rasterio.open(in_situ_val_tif, 'w+', **out_meta)
dst_arr = dst.read(1)
# This is where we create a generator of geom, value pairs to use in rasterizing
geom_col = in_situ_gdf.geometry
code_col = in_situ_gdf[field_classif_code].astype(int)
shapes = ((geom,value) for geom, value in zip(geom_col, code_col))
in_situ_arr = features.rasterize(shapes=shapes,
fill=no_data,
out=dst_arr,
transform=dst.transform)
dst.write_band(1, in_situ_arr)
print(f'Rasterize is done : {in_situ_val_tif}')
# Close rasterio objects
src.close()
dst.close()
else:
print('CRS are different --> repoject in-situ data shapefile with "to_crs"')
Raster template file : /export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/CLASSIF/NAMUR_2020_classif_RF_feat_2_grp_A.tif The CRS of in situ data is : 32631 The CRS of raster template is : 32631 CRS are the same Rasterize starts : /export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/IN_SITU/NAMUR_2020_IN_SITU_ROI_VAL.shp Rasterize is done : /export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/IN_SITU/NAMUR_2020_IN_SITU_ROI_VAL.tif
y_pred and y_true¶# Open in-situ used for validation
src = rasterio.open(in_situ_val_tif, "r")
val_arr = src.read(1)
src.close()
# Open classification map
src = rasterio.open(classif_tif, "r")
classif_arr = src.read(1)
src.close()
# Get the postion of validation pixels
idx = np.where(val_arr == no_data, 0, 1).astype(bool)
# Ground truth (correct) target values
y_true = val_arr[idx]
# Estimated targets as returned by a classifier.
y_pred = classif_arr[idx]
Sometimes, some classes do not appear in the classification map, they are not predicted by the Random Forest.
This means that some classes in y_true don't appear in y_pred.
classes_all = sorted(np.unique(y_true))
classes_pred = sorted(np.unique(y_pred))
classes_missing = set(y_true) - set(y_pred)
print(f'{len(classes_missing)} classes are missing in the classification (y_pred) : {classes_missing} \n')
print(f'All training classes :\n {classes_all}')
print(f'All predicted classes (at least once) :\n {classes_pred}')
0 classes are missing in the classification (y_pred) : set() All training classes : [3, 6, 8, 9, 17, 21, 22, 111, 112, 115, 117, 119, 121, 143, 151, 181, 192] All predicted classes (at least once) : [3, 6, 8, 9, 17, 21, 22, 111, 112, 115, 117, 119, 121, 143, 151, 181, 192]
lut_df = pd.read_excel(s4s_lut_xlsx)
classes_name = lut_df[lut_df[field_classif_code].isin(classes_all)].sort_values(field_classif_code)[field_classif_name].drop_duplicates().to_list()
for code,name in zip(classes_all, classes_name):
print(f'{code} - {name}')
3 - Grassland and meadows 6 - Forest 8 - Build-up surface 9 - Water bodies 17 - Leguminous crops 21 - Fruits trees 22 - Vineyards 111 - Wheat 112 - Maize 115 - Barley 117 - Oats 119 - Other cereals 121 - Leafy or stem vegetables 143 - Other oilseed crops 151 - Potatoes 181 - Sugar beet 192 - Fibre crops
cm = confusion_matrix(y_true, y_pred)
cm_df = pd.DataFrame(cm)
cm_values = cm_df.values
cm_df.columns = classes_all
cm_df.index = classes_all
cm_df.to_csv(cm_csv, index=True, sep=',')
display(cm_df)
| 3 | 6 | 8 | 9 | 17 | 21 | 22 | 111 | 112 | 115 | 117 | 119 | 121 | 143 | 151 | 181 | 192 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 3 | 17445 | 298 | 20 | 0 | 5 | 56 | 20 | 62 | 0 | 1 | 3 | 8 | 0 | 0 | 0 | 0 | 0 |
| 6 | 37 | 746 | 30 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 8 | 52 | 14 | 422 | 13 | 0 | 3 | 0 | 20 | 0 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 9 | 0 | 0 | 26 | 30 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 17 | 151 | 0 | 7 | 0 | 1116 | 0 | 0 | 20 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 8 |
| 21 | 686 | 0 | 0 | 0 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 22 | 14 | 0 | 0 | 0 | 0 | 40 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 111 | 90 | 3 | 47 | 0 | 14 | 0 | 4 | 20100 | 0 | 1059 | 65 | 750 | 0 | 2 | 0 | 1 | 95 |
| 112 | 28 | 0 | 25 | 0 | 2 | 0 | 0 | 1 | 6344 | 0 | 3 | 10 | 4 | 7 | 87 | 22 | 2 |
| 115 | 199 | 0 | 4 | 0 | 0 | 0 | 0 | 1509 | 0 | 2726 | 136 | 401 | 0 | 21 | 0 | 0 | 0 |
| 117 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 11 | 0 | 0 | 88 | 0 | 0 | 0 | 0 | 0 | 142 |
| 119 | 12 | 0 | 11 | 0 | 0 | 0 | 0 | 4005 | 0 | 0 | 112 | 455 | 0 | 0 | 0 | 0 | 45 |
| 121 | 5 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 738 | 0 | 0 | 0 | 23 | 0 | 0 | 19 | 0 |
| 143 | 51 | 0 | 0 | 0 | 0 | 0 | 0 | 18 | 0 | 24 | 0 | 0 | 0 | 557 | 0 | 0 | 0 |
| 151 | 0 | 1 | 6 | 0 | 0 | 0 | 1 | 0 | 271 | 0 | 0 | 4 | 6 | 0 | 3538 | 1035 | 0 |
| 181 | 2 | 14 | 1 | 0 | 4 | 2 | 1 | 0 | 389 | 0 | 0 | 0 | 166 | 0 | 21 | 4950 | 0 |
| 192 | 0 | 0 | 2 | 0 | 0 | 0 | 0 | 605 | 0 | 0 | 40 | 45 | 0 | 0 | 0 | 0 | 59 |
# invert z idx values
z = cm[::-1]
x = classes_name
y = x[::-1].copy() # invert idx values of x
# change each element of z to type string for annotations
z_text = [[str(y) for y in x] for x in z]
# set up figure
fig = ff.create_annotated_heatmap(z,
x=x,
y=y,
annotation_text=z_text,
colorscale='spectral',
reversescale=True)
# add title
fig.update_layout(title_text=f"Confusion Matrix - {site}, {year}")
# adjust margins to make room for yaxis title
#fig.update_layout(margin=dict(t=200, l=200))
#fig.update_xaxes(tickfont_size=20)
#fig.update_yaxes(tickfont_size=20)
#fig.update_layout(font_size=25)
# add colorbar
#fig['data'][0]['showscale'] = True
fig.show()
fig.write_html(cm_html, full_html=False)
If you decide that you are not interested in the scores of classes that were not predicted, then you can explicitly specify the classes you are interested in (which are labels that were predicted at least once).
acc_metrics_str = classification_report(y_true,
y_pred,
target_names=classes_name,
labels=classes_all,
digits=3)
print(acc_metrics_str)
precision recall f1-score support
Grassland and meadows 0.929 0.974 0.951 17918
Forest 0.693 0.918 0.790 813
Build-up surface 0.702 0.802 0.749 526
Water bodies 0.698 0.536 0.606 56
Leguminous crops 0.976 0.856 0.913 1303
Fruits trees 0.000 0.000 0.000 688
Vineyards 0.000 0.000 0.000 54
Wheat 0.763 0.904 0.827 22230
Maize 0.819 0.971 0.889 6535
Barley 0.715 0.546 0.619 4996
Oats 0.197 0.364 0.255 242
Other cereals 0.272 0.098 0.144 4640
Leafy or stem vegetables 0.116 0.029 0.047 785
Other oilseed crops 0.949 0.857 0.901 650
Potatoes 0.970 0.728 0.832 4862
Sugar beet 0.821 0.892 0.855 5550
Fibre crops 0.168 0.079 0.107 751
accuracy 0.807 72599
macro avg 0.576 0.562 0.558 72599
weighted avg 0.774 0.807 0.783 72599
oa = accuracy_score(y_true, y_pred)
oa = round(oa*100, 2)
print(f'Overall Accuracy : {oa}%')
Overall Accuracy : 80.72%
acc_metrics_dict = classification_report(y_true, y_pred,target_names=classes_name, output_dict=True)
am_df = pd.DataFrame.from_dict(acc_metrics_dict).round(3)
am_df = am_df.iloc[:,:-3]
#am_df = pd.concat([am_df, nb_df])
#am_df = am_df.sort_values(by='pix_count', ascending=False, axis=1)
class_name = am_df.columns.to_list()
precision = am_df.loc['precision'].to_list()
recall = am_df.loc['recall'].to_list()
f1_score = am_df.loc['f1-score'].to_list()
fig = go.Figure(data=[
go.Bar(name='Precision', x=class_name, y=precision, text=precision, textposition='auto'),
go.Bar(name='Recall', x=class_name, y=recall, text=recall, textposition='auto'),
go.Bar(name='F1-score', x=class_name, y=f1_score, text=f1_score, textposition='auto')
])
# Change the bar mode
fig.update_layout(title_text=f'Accuracy Metrics - {site}, {year}',
barmode='group')
#fig.update_xaxes(tickfont_size=30)
fig.update_yaxes(tickfont_size=10, range=[0,1])
#fig.update_layout(xaxis_title=None, font_size=10)
#fig.update_layout(legend=dict(font=dict(size=25)))
fig.show()
fig.write_html(am_html, full_html=False)